Strain-induced kinetics of intergrain defects as the mechanism of slow dynamics in the 
nonlinear resonant response of humid sandstone bars 
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A closed-form description is proposed to explain nonlinear and slow dynamics effects exhibited by 
sandstone bars in longitudinal resonance experiments. Along with the fast subsystem of longitudi- 
nal nonlinear displacements we examine the strain-dependent slow subsystem of broken intergrain 
and interlamina cohesive bonds. We show that even the simplest but phenomenologically correct 
modelling of their mutual feedback elucidates the main experimental findings typical for forced lon- 
gitudinal oscillations of sandstone bars, namely, (i) hysteretic behavior of a resonance curve on both 
its up- and down-slopes, (ii) linear softening of resonant frequency with increase of driving level, 
and (iii) gradual recovery (increase) of resonant frequency at low dynamical strains after the sample 
was conditioned by high strains. In order to reproduce the highly nonlinear elastic features of sand- 
stone grained structure a realistic non-perturbative form of strain potential energy was adopted. 
In our theory slow dynamics associated with the experimentally observed memory of peak strain 
history is attributed to strain-induced kinetic changes in concentration of ruptured inter-grain and 
inter-lamina cohesive bonds causing a net hysteretic effect on the elastic Young's modulus. Finally, 
we explain how enhancement of hysteretic phenomena originates from an increase in equilibrium 
concentration of ruptured cohesive bonds that are due to water saturation. 
PACS numbers: 05.45.-a, 62.40.+i, 83.80.Fg, 46.05. +b 



Apart from their excellent static characteristics as 
building materials, sandstones have been shown to 
demonstrate a number of unexpected and even surprising 
dynamical properties [1-5] . Here we consider the numer- 
ous experimental results on nonlinear resonant response 
exhibited by sandstone rods in forced longitudinal oscil- 
lations that appear even at extremely small forcing lev- 
els and consequently at small dynamic strains [1-5] . The 
most intriguing nonlinear feature is slow dynamics, which 
is defined here as long-term (minutes to hours) change of 
elastic properties in response to a disturbance such as 
dynamic and static strain or temperature. 

Specifically, we have to underline that in the vicin- 
ity of bar resonant frequency the longitudinal alternat- 
ing drive produces strong essentially non-trivial nonlinear 
responses: 1) At high drive levels the effective width of 
resonance curves depends on the direction of frequency 
sweep; it is narrower for upward sweeps (i.e., from lower 
to higher frequencies) than at downward sweeps (i.e., 
from higher to lower frequencies) [1-5] . This effect proves 
to be a typical manifestation of slow dynamics and can 
be treated as hysteresis both on low- and high-frequency 
slopes of a resonance curve. 2) The resonance peak is 
shifted toward lower frequency almost linearly with in- 
crease of driving amplitude [1, 4]. 3) Other evidence 
of slow dynamics comprises gradual recovery (increase) 
of resonant frequency to its original value as defined at 
extremely low drive level after the sample has been con- 
ditioned at high drive level [3,5]. 

These facts cannot be understood in the framework of 
standard theories of resonant nonlinear response [6] and 
imply memory of peak strain history [2] . Some aspects of 
the problem have been explained by the interpretation of 



Guyer, McCall and Van Den Abeele [7] in the framework 
of a McCall-Guyer quasistatic model [8]. This approach 
uses the concept of auxiliary hysteretic elements that al- 
low the introduction of an additional nontrivial nonlinear 
term into the dynamical equation for the field of longi- 
tudinal displacements. However, this theoretical treat- 
ment lacks completeness in that it initially neglects the 
dynamics of hysteretic elements and postulates tempo- 
ral evolution of amplitude-frequency characteristic (the 
key point of claimed results) to be developed afterwards. 
Although Capogrosso-Sansone and Guyer recently sug- 
gested a dynamical realization of the McCall-Guyer qua- 
sistatic model [9], evaluating its adequacy to explain ex- 
perimental data turns out to be difficult. 

In this communication we omit the idea of auxiliary 
hysteretic elements as the sole approach for treating all 
peculiar hysteretic phenomena and call attention to an al- 
ternative notion used by Davydov and Ermakov for the 
description of bistability in nonlinear resonant tunnel- 
ing of electrons through a set of potential barriers [10]. 
Their approach consists of explicit but physically moti- 
vated separation of given physical system into two non- 
linear subsystems, namely fast and slow subsystems with 
mutual coupling taken into account. 

For sandstone bars we identify the fast subsystem with 
the field of rapid longitudinal displacements while the 
slow subsystem represents the concentration of defects in 
intergrain contact bonds. In doing this we bear in mind 
that, because of preferable vapor condensation onto sur- 
faces with greater concave curvature [11], the sandstone 
pore structure [4, 11] retains some residual pore water 
[11], and its impact on the resonant properties of rock is 
crucial [12, 13]. Thus, thermodynamical estimations ap- 
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plied to porous rocks show that intergrain cohesive forces 
become weaker in the presence of water [14], that agrees 
with an alternative conception of swelling pressure [13, 
15]. This treatment is supported by recent experiments 
[13] establishing an abrupt decrease in Young's modulus 
within the first twenty percent interval of water satu- 
ration (i.e., until the degree that void surfaces become 
completely wet). We could additionally invoke ordinary 
capillary forces [13] or hydrolysis of silicon-oxygen-silicon 
bond-chains [16] in our consideration. However, either of 
these mechanisms also leads to softening of Young's mod- 
ulus with saturation increase. Here the significant issue 
is apparently not in excessive (presumably unclaimed) 
detailing all plausible mechanisms that might modify the 
Young's modulus in a qualitatively similar way, but in 
their reasonable concise formalization by means of a min- 
imal number of slow fields. 

According to Kosevich [17] the equilibrium concentra- 
tion of defects associated with a stress a is given by the 
formula 



c a = cq exp (va/kT) , 



(1) 



where k and T are the Boltzmann constant and the abso- 
lute temperature, respectively, and the parameter v > 
stands for a typical volume accounting for a single defect 
and characterizes the intensity of dilatation. The equi- 
librium concentration of defects in an unstrained bar cq 
has to be some function of both temperature T and water 
saturation s. In order to describe strain-induced changes 
in nonequilibrium concentrations of defects c, we assume 
that at any instant of time t the concentration c must 
evolve to its would-be equilibrium value c CT , where the 
stress u in (1) is applied at the same instant. Supposing 
the distributions of activation barriers for defect annihila- 
tion U and activation barriers for defect creation W to be 
uniform respectively over the ranges Uq < U < Uq + U+ 
and Wo < W < W + W+ with U , U+ and W , W+ be- 
ing insensitive to the choice of bar cross-section, we will 
deal with the density of defect concentration g governed 
by the following kinetic equation 

dg/dt = - [ M % - g a ) + v6(g a - g)] (g - g a ). (2) 

Here fi = fi cxp(— U/kT) and v = v a cxp(— W/kT) are 
the rates of defect annihilation and defect creation re- 
spectively, g a = c a /U + W + , and 9(z) designates the 
Heaviside step-function. The quantities g and c arc re- 
lated by the simple definition 



rU +U+ rW +W+ 

f= / dU dW-g. 



(3) 



Under the tensile load there is an immense number of 
spatial ways for the intergrain cementation contact to be 
cleaved with the same basic result: the creation of crack. 
The similar scenario is true also for the already existed 
balanced crack to be further expanded. On the contrary 
under the compressive load the crack ones emerged has 



only one spatial way to be annihilated or contracted. 
These observations are the principal ones and imply the 
huge disparity v 3> p between the priming rates vq and 

notwithstanding the generic cohesive properties of ce- 
mentation material. Moreover, because of possible frag- 
mentation of cementation material and/or water interca- 
lation between the opposite faces of crack we can expect 
the typical value of barrier U to exceed that of barrier W. 
In combination all these factors might sustain predomi- 
nantly even the more immense disparity v 3> fj, between 
the actual rates v and fi of defect creation and defect 
annihilation, apparently comprising many orders, and as 
a result provide the physical mechanism that breaks the 
symmetry of system response to an alternating external 
drive and acts as a sort of soft ratchet or leaky diode. 

To express the evolution equation 



cPu 



da d 

dx dx 



dT 



d(d 2 u/dxdt) 



(4) 



for the field of longitudinal displacements u we choose 
the stress-strain relation in the form 



a = 



Esechrj 



(r — a)[cosh?y du/dx + l] a+1 

-Esechf? 
(r — a) [cosh rj du/dx + l] r+1 



(5) 



which at r > a > allows one to block the bar compress- 
ibility at strains du/dx tending to +0 — scchry. To put it 
differently the parameter sechr; is reserved for the typical 
thickness of intergrain cementation contact divided by 
the typical distance between the centers of neighbouring 
grains. The dissipative function T we take in the form 



T= (7/2) [d 2 u/dxdt\ 



(6) 



giving rise to Stokes internal friction [18]. Here x denotes 
the longitudinal Lagrange coordinate of the bar sample. 
The quantities p and 7 are respectively the mean den- 
sity of sandstone and the coefficient of internal friction 
in an elastic subsystem. We ignore their dependences 
on temperature and water saturation assuming that the 
main effect is manifested through the linear decrease of 
Young's modulus E with the concentration of defects 



E= (1-c/c 



<E, 



(7) 



Here c cr and E + are the critical concentration of defects 
and the maximal possible value of Young's modulus, re- 
spectively. Both of these parameters we also take to be 
independent of temperature and water saturation. 

Typical resonant response experiments [1-5] corre- 
spond to the kinematic excitation [19] of a bar sample, 
which we associate with the following boundary condi- 
tions 



u(x = 0\t) = D(t) cos + cJtw(t)^ (8) 

(9) 



du , , . . 
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where L is the sample length and t > 0. The driving 
amplitude D(t) is assumed to be basically constant ex- 
cept for the moments when the driving system is switched 
on, is switched into another constant driving level, or is 
switched off. The time dependence of cyclic driving fre- 
quency u(t) in turn is prescribed by the chosen regime of 
frequency sweep. Initial conditions are given in the form 



u{x\t = 0) = 
g(x\t = 0) = c Q /U+W+ 



du , 



0) = (0 < x < L) (10) 
(0 < x < L). (11) 



When experimental data for Young's modulus in un- 
strained samples are obtainable from the resonant re- 
sponse experiments by the use of a low amplitude proto- 
col (driving amplitude D very small and negligible strain- 
induced defects), we can compare them with values taken 
from equation (7) at c = cq in order to fit the equilib- 
rium concentration of defects c as a function of T and 
s by some extrapolation formula. In particular, relying 
upon the Sutherland temperature extrapolation [20] and 
analyzing temperature dependent data at zero saturation 
[21] and saturation dependent data at room temperature 
[13] for Berea sandstone, we suggest the formula 



Co = c c 



T 

T7r 



cosh 2 a — exp ( — S \ sinh 2 a 



(12) 



with the following fitting parameters T cr =1475 K, 
cosh 2 a = 16, /3 = 10. Here the saturation can vary 
within the interval < s < 1. Ats^O this approx- 
imation is expected to work at least for temperatures 
exceeding the freezing-point of pore water. 

Computer modelling of nonlinear and slow dynamics 
effects was performed in the vicinity of the resonant fre- 
quency /o(2), which we understand as the second fre- 
quency (I = 2) in the fundamental set 



/o(0 = 



21 - 1 




(1 = 1,2,3,...) (13) 



given by the linear theory of kinematic excitation for zero 
dissipation 7 = 0. 

Figure 1 shows typical resonance curves, i.e. depen- 
dences of response amplitude R (calculated at x = L) on 
drive frequency / = oj/2tt, at successively higher drive 
amplitudes D. The continuous lines correspond to the 
conditioned resonance curves calculated after two fre- 
quency sweeps were performed at each driving level in 
order to achieve repeatable hysteretic curves. The dot- 
ted line illustrates an unconditioned curve obtained with- 
out any preliminary conditioning. Arrows on the three 
highest curves indicate sweep directions. For the sake of 
definiteness the results of the computer simulation were 
adapted to the experimental conditions for the data ob- 
tained by Ten Cate and Shankland for Berea sandstone 
[2] . In particular, the ratio E + / p was estimated by means 



of relationships (13) and (12) with the second order fre- 
quency, bar length, temperature and saturation as fol- 
lows: / (2)=3920Hz, L = 0.3m, T = 297 K and s = 0.25. 
The ratio 7/p characterizing internal friction was chosen 
from the best fit to low amplitude theoretical (Figure 
1) and experimental [2] resonance curves using the qual- 
ity factor Q from resonance lineshape. The parameters 
(i exp(-U /kT) = Is- 1 and U+/k = 2525 K determin- 
ing the character of slow relaxation were estimated ac- 
cording to the experimental measurements of decay of 
acceleration at fixed frequency [2] and the observations 
of recovering resonant frequency as a function of time [5] . 
The combination of parameters vE + /kcoshr) = 275 K 
was chosen to quantitatively reproduce the hysteretic 
phenomena in the sweep regimes typical for the actual ex- 
periments [2] . The nonlinearity parameter cosh 77 = 2300 
was estimated to map the true asymmetry of experimen- 
tal resonance curves [2]. Other parameters appearing in 
the stress-strain relation (5) have been adopted as fol- 
lows: r = 4, a = 2. 




Driving frequency (Hz) 

FIG. 1. Resonance curves j = 0,1,2,3,4,5 at suc- 
cessively higher driving amplitudes Dj : Dj/L = 
= 3.8[?'(1 — Sjo) + 0.5(5jo] ■ 1CP 8 . Continuous lines denote con- 
ditioned curves; the dotted line represents the unconditioned 
curve. Arrows on the highest curves indicate sweep direc- 
tion. The time to sweep back and forth within the frequency 
interval 3700Hz^4100Hz is chosen to be 120 s. 

We would like to stress that through the drop of equi- 
librium concentration cq our theory is able to catch the 
drastic suppression of hysteresis with decrease of water 
saturation (see Eq. (12)). This conclusion has been con- 
firmed by direct computation (not shown). Simultane- 
ously we have observed a monotonic increase in quality 
factor Q with saturation decrease, i.e., precisely the well 
documented tendency in experiments [12]. 

Figure 2 compares the shifts of resonant frequency as 
functions of driving amplitude at two essentially different 
values of dilatation parameter v while other parameters 
were kept the same as for Figure 1. Thus curve 1 calcu- 
lated at vE + /k cosh -q = 275 K, when the strain-induced 
feedback between the slow and the fast subsystems is 
substantial, demonstrates the almost linear dependence 
typical for materials with nonclassical nonlinear response, 



3 



i.e., materials possessing all the basic features of slow dy- 
namics. On the other hand, curve 2 calculated at v = 0, 
when the strain-induced excitation of slow subsystem is 
absent and hence the mutual feedback between the slow 
and the fast subsystems is totally broken, demonstrates 
the almost quadratic dependence typical for the materi- 
als with classical nonlinear response. 
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FIG. 2. The shift f r — fo of resonant frequency f r from 
its asymptotic value fo as the function of driving amplitude 
D for the hysteretic nonlinear material (curve 1) and for the 
classical nonlinear material with v = (curve 2). 



of excitation and (ii) the shift caused by the effect of the 
slow subsystem. However, only the second part might ac- 
tually be registered during the recovery process, because 
the first one vanishes almost instantaneously when the 
conditioning was switched off. Hence, the whole charac- 
ter of recovery should inevitably be governed by the slow 
kinetics responsible for restoration of intergrain bonds 
[5]. From Figure 3 we clearly see that the suggested ki- 
netic equation for the density of defect concentration (2) 
supplemented by the simple definition of total concentra- 
tion (3) and the reasonable relationship between Young's 
modulus and a concentration of defects (7) yields the 
very wide time interval 10 < (t — t c )/t < 1000 of log- 
arithmic recovery of the resonant frequency in complete 
agreement with experimental results [5]. Here t c is the 
moment when conditioning switches off and to=l s. 

This work was carried out within the framework of 
project No 1747 supported by the STCU. O.O.V. ac- 
knowledges support from the National Academy of Sci- 
ences of Ukraine (Grant 0102U002332). We would like 
to thank Paul A. Johnson for his unflagging interest to 
the process of multi-stage amendments of the model. 




FIG. 3. Time-dependent recovery of resonant frequency f r 
to its asymptotic value /o after the large conditioning drive 
has been removed. Curves 1, 2, 3 correspond to successively 
higher saturations s — 0.05, s = 0.15, s — 0.25. The fre- 
quency shift f r — fo is normalized by both the asymptotic fre- 
quency fo and the unitless response amplitude R/L attained 
at conditioning resonance. 

Finally, Figure 3 shows the gradual recovery of reso- 
nant frequency f r to its maximum limiting value fo after 
the bar has been subjected to high amplitude condition- 
ing and conditioning is stopped. In the computer simula- 
tion we have plotted three different curves corresponding 
to three different saturations with all other model param- 
eters adopted earlier for Figure 1 being preserved. The 
total shift of resonant frequency f r — fo consists of two 
physically different parts, namely (i) the traditional dy- 
namic shift caused by strain nonlinearity at high levels 
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